%This script generate the computational simulation grid 
%for each ring-array location.
clear all;
%parameters
dx = 0.2;  %[mm] pixel size 
PML_size = 20;  %  
buffer_layer_size = 20; %the size of apodization layers, as described in paper Appendix B
z_interval = 1.8; %[mm] the scanning interval of two adjacent ring-array
z_voxel_shift = z_interval/dx;

type = 'B'; %select phantom type B or C


if type == 'B'
    prelix = 'scatter';
elseif type == 'C'
    prelix = 'hetero';
else
    error('The selected phantom type must be B or C') 
end

% load phantom
phantom_name = ['./',type,'/phantom.mat'];
load(phantom_name);
assert(all(size(sos) == [1280,1280,194]),'check the input phantom size')
assert(all(size(aa) == [1280,1280,194]),'check the input phantom size')
assert(all(size(density) == [1280,1280,194]),'check the input phantom size')

Nz = size(sos,3)-z_voxel_shift*2;
alpha = get_alpha(buffer_layer_size, Nz); % get the apodization weights
c_bg = sos(1); d_bg=density(1); a_bg = aa(1); % get background value;


for ring = [-1,0,1]
    % extract thin-slab for certain ring-array
    sos_temp = sos(:,:,(ring+1)*z_voxel_shift+1:end+(ring-1)*z_voxel_shift);
    sos_temp = alpha.*c_bg+sos_temp.*(1-alpha); % smooth the top and bottom boundary
    density_temp = density(:,:,(ring+1)*z_voxel_shift+1:end+(ring-1)*z_voxel_shift);
    density_temp = alpha.*d_bg+density_temp.*(1-alpha); %smooth the top and bottom boundary
    aa_temp = aa(:,:,(ring+1)*z_voxel_shift+1:end+(ring-1)*z_voxel_shift);
    %extend the media for PML layers
    sos_slab = c_bg*ones(1280,1280,Nz+PML_size*2);
    density_slab = d_bg*ones(1280,1280,Nz+PML_size*2);
    aa_slab = a_bg*ones(1280,1280,Nz+PML_size*2);
    sos_slab(:,:,PML_size+1:end-PML_size)=sos_temp;
    density_slab(:,:,PML_size+1:end-PML_size)=density_temp;
    aa_slab(:,:,PML_size+1:end-PML_size)=aa_temp;
    out_name = ['./',type,'/phantom_ring_{',num2str(ring),'}.mat'];
    save(out_name, 'sos_slab','density_slab','aa_slab','-v7.3');
end

function alpha=get_alpha(buffer_size,Nz)
% function to compute the apodization weights
    alpha =zeros(1,1,Nz);
    alpha(1,1,1:buffer_size) = cos(pi/2*(1:buffer_size)/buffer_size).^2;
    alpha(1,1,Nz-buffer_size+1:Nz) = cos(pi/2*(buffer_size:-1:1)/buffer_size).^2;
end